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Abstract 

We present a data mining technique for the analysis of multichannel oscillatory 
timeseries data and show an application using poloidal arrays of magnetic sensors 
installed in the H-1 heliac. The procedure is highly automated, and scales well to 
large datasets. The timeseries data is split into short time segments to provide 
time resolution, and each segment is represented by a singular value decomposition 
(SVD). By comparing power spectra of the temporal singular vectors, singular values 
are grouped into subsets which define fluctuation structures. Thresholds for the 
normalised energy of the fluctuation structure and the normalised entropy of the 
SVD are used to filter the dataset. We assume that distinct classes of fluctuations 
are localised in the space of phase differences A^/^(n, n + 1) between each pair of 
nearest neighbour channels. An expectation maximisation clustering algorithm is 
used to locate the distinct classes of fluctuations, and a cluster tree mapping is used 
to visualise the results. 
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1 Introduction 



The motivation for the present work arose from the analysis of fluctuations 
in magnetically confined plasma during parameter scans in the H-1 flexible 
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heliac [T][2] . The H-1 heliac is a three field-period hehcal axis stellarator [3] 
with major radius R = Im, minor radius (r) = 0.2m and a finely tunable 
magnetic geometry. 

Experimental scans through plasma configurations via the geometric param- 
eter Kh, which controls the rotational transform t (twist of the magnetic field 
lines) and shear t'(radial derivative of rotational transform), have produced 
diverse spectra of magnetohydrodynamic (MHD) activity. The MHD activity 
is recorded via two toroidally separated poloidal arrays of Mirnov coils (induc- 
tion solenoids) which sample dB/dt locally. In the example dataset presented 
here, 28 Mirnov coils are used for 92 distinct plasma configurations, resulting 
in more than 100,000 short time Fourier spectra. 

The data mining process used to reduce this dataset is described in the follow- 
ing sections. In section [2] we explain the preprocessing stage, which includes 
filtering and mapping into a high dimensional phase space. In section [3] the 
clustering algorithm for distinguishing classes of fiuctuations is described, fol- 
lowed by a demonstration of a visualisation procedure. A discussion of some 
important aspects of the procedure follows in section HI 



2 Preprocessing 

2.1 Data preparation 

We assume that each set of timeseries data can be represented as a A^c x Ng 
matrix: 



V SnM SnA^o + SN,{to + 2r) ... SN,{to + N,t) J 

where r is the inverse of the sampling frequency, Nc is the number of channels 
and Ns is the number of samples. In our example dataset, the signal amplitudes 
depend on the plasma-coil distance which is a function of the plasma shape 
(magnetic configuration) controlled by Kh- To reduce any configurational bias 
on S we normalise each channel to its variance. 

To achieve time resolution At, we split S into short time segments S with 
shape Nc x iV^, where N'^ = At/r. We also assume there are an arbitrary 
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number of S relating to the same system, e.g.: an experiment repeated under 
different conditions. At this stage there is no need to distinguish between the 
S from different S, although we implicitly retain sufficient information to map 
the S back to their original parameter sets. 



2.2 The singular value decomposition 

Each S is represented by a singular value decomposition (SVD) [1] 

S = UAV* (2) 



where the columns of U and V contain the spatial (topo) and temporal 
(chrono) singular vectors respectively, V* denotes the conjugate transpose of 
V , and the diagonal elements of A are the Na = min(A'^c, N'^) non-negative sin- 
gular values. The set of topos (chronos) are an orthonormal basis of R^='-^«^ 
The convention is for the singular values to be sorted in decreasing monotonic 
order meaning that A is independent of the ordering of the channels within 
5*. Shown in figure [1] are singular values from a typical HI Mirnov dataset. 
From the chronos power spectra we see that there are two dominant modes, 
each with two singular values suggesting they are both travelling waves, as 
discussed below. We also see that the variance-normalisation of each channel 
degrades the signal to noise ratio of the system, which can also be described 
in terms of the normalised entropy. 

We calculate the normalised entropy H of the singular values in A: 



where pk is the dimensionless energy: 

P. = |, E = Y.al (4) 

fc=i 



The low entropy case {H —>■ 0) occurs when the system is well ordered. To 
some extent the scalar quantity H can be used as a measure of how physically 
interesting the signals in S are without any further investigation into the 
structure of S, though care must be taken with this interpretation. A standing 
wave in a system with no noise has only one non-zero singular value ag = 1 
giving H = 0, whereas a travelling wave requires two singular values so H > 0. 
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Fig. 1. Example of chronos power spectra and singular values. Singular values from 
both normalised (o) and unnormalised (x) S are shown. C0,C1,...,C5 denote 
the chronos from the normalised singular value 0, 1, . . . , 5. There are two distinct 
modes, one at / ~ 45 kHz described by SVO and SVl; the other, weaker, signal is 
at / ~ 29 kHz and is described by SV2 and SV3. The data is from HI shot #58122 
at 31 < t < 32 ms. 

2.3 Fluctuation structures 



Recognising that a travelling wave structure consists of a pair of singular values 
which naturally belong together we group similar singular values, defining a 
fluctuation structure a as a subset of singular values which have chronos with 
similar power spectra. We measure the similarity between two chronos Ci and 
C2 with the normalised average of the cross-power spectrum 7ci,c2- 

^ g(ci,C2)^ 

^^^•^^ G{c^,c,)G{c„c,y 

where G{a, h) = (|jF(a)jF*(6)|), jFisthe Fourier transform, and (. . .) represents 
the spectral average. 

When allocating singular values to fluctuation structures, the observation: 

and 7a,c > 7mm 76, c ^ 7™n) (6) 



suggests that we should not simply seek to require 7a,b > 7mm for each pair 
of singular values a, b within a structure, instead we follow the process in 
algorithm 1. In so doing we therefore require that each constituent singular 
vector has sufficient 7 with the dominant singular vector of the structure. 

Various possible fluctuation structures for the dataset of figure [1] are shown 
in figure [2] as a function of the threshold value 7mm- At 7mm = 0, all singular 
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while Number of unallocated singular values > do 

Define a new fluctuation structure as an empty set of singular values: 
ai = {} 

Denote the largest unallocated singular value by 
for Every unallocated singular value do 
if 7^,g > 'jmin then 

Allocate to fluctuation structure at 
end if 
end for 
end while 

Algorithm 1: Building fluctuation structures from singular values aj. The 
largest unallocated singular value will always be allocated to because 
7?,C = 1- 
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Fig. 2. The possible fluctuation structure groupings according to their energies as 
defined by algorithm 1 through tlie range of 7mm- The dataset is the same as in figure 
m We see that 7mm ^ 0.4 aUows unrelated singular values to be included within a 
fluctuation structure, whereas with 7mm > 0.87 algorithm 1 will not recognise the 
similarity between a2 and 03. 



values are grouped together as a single fluctuation structure, while at 7mm = 1 
each fluctuation structure contains one singular value. The key features are 
the two fluctuation structures = Wo, 0.1} and ai = {02,03} which coexist 
for 0.50 < 7mm < 0.87. After application of such analysis to a suitably sized 
sample of short time segments, a threshold of 7mm = 0.7 was found to be 
appropriate for our dataset. 
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2.J^ Data filtering 



Filters are applied to the dataset in order to reduce its size and to remove 
noise. Two values which can be used to quantify the quality of the data are 
the normalised energy p and normalised entropy H. The normalised energy p 
of a fluctuation structure is defined as the sum of the normalised energies of 
its constituent singular values from equation HI The nature of these thresholds 
is quite different; H thresholds will act upon the entire short time segments, 
whereas p thresholds affect individual fluctuation structures. 

Using a normalised energy threshold value p' allows filtering out of low energy 
noise. As seen in our example dataset (figure [3]), there appears to be a clear 
distinction between higher energy fluctuations {p > 0.6) and lower energy 
noise (p < 0.2). The use of H thresholds is not always appropriate, especially 
if spectra are present with several distinct fluctuations; in such cases if the 
dataset needs to be reduced in size it is preferable to simply use a random 
subset of the data. Entropy filtering is generally more useful when a significant 
number of S contain only noise. 

An alternative to using a hard p threshold uses an energy threshold defined as 
a fraction of the possible range of normalised energy for a given singular value, 
before constructing fluctuation structures. This method is more sensitive to 
modes which have reduced pk due to the coexistence of other modes. The 
n**^ largest singular value in a given short time segment has, by definition, a 
maximal normalised energy of 1/n. We then apply a factor p* , where < p* < 
1, and, starting with the smallest singular value, retain singular values with 
Pn > p*/n. Any point retained brings in all larger singular values from the same 
time segment, trumping the Pn > p*/n condition for lower n. The requirement 
for bringing in larger singular values is easily understood by considering the 
case of a mode having two singular values of almost equal energy, it is possible 
for the lower energy value to exceed the threshold with the higher value below 
the threshold. 



2.5 Mapping of fluctuation structures into Aip-space 

We regard each fluctuation structure as a point in the space [— vr, ir]'^", an A^^c- 
dimensional torus of length 2tt which we will call Aip-space. In this application 
ip represents the electrical phase of the reconstructed fluctuation structure 
at the positions of the coils. Fluctuation structures which are close in A^- 
space can be considered the same type. This interpretation arises from the 
expectation that the possible waves have various phase velocities and mode 
numbers due to the periodic boundary conditions of the physical plasma torus 
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Fig. 3. A 10% random sample of a datasct. The left panel shows p and H for the 
fluctuation structures. The middle panel shows the number of fluctuation structures 
Na within Sp = 0.01. The right panel shows pN^, which is effectively the density 
of normalised energy; while this is not physically meaningful because the normali- 
sation factor is dependent on short time segment, it is a useful guide to the energy 
distribution among fluctuation structures. 

in which they propagate. It is also applicable to the more general case where 
waves arc spatially localised within the system and do not have well defined 
mode numbers. 

For each fluctuation structure ai we take the inverse SVD to get Sf. 

Si = UAiV*, (7) 



where the elements of A not in ai are set to zero to form Ai. The rows of the 
matrix Si contain the timeseries relating to ai for each channel. In general, 
the power spectra of the topos in ai are peaked around a single frequency uji. 
The phase differences Aipajii^^ = ooi) between channels a and h evaluated at 
uj = uji are used to define the coordinates in A-i/'-space. Using phase differences 
between each pair of channels would result in a |iVc(iVc — l)-dimensional space; 
instead we use the A^c "dimensional space of only nearest neighbour channels. 
Note that in our example dataset the actual phase difference between channels 
depends on so we map the phase differences to a coordinate system which 
is independent of k/j, namely the K/j-averaged magnetic angles of the Mirnov 
coils. 
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Fig. 4. The preprocessed dataset: In the upper panel, fluctuation structures are 
mapped to frequency and the magnetic geometry parameter k/j, marker size and 
colour are proportion to the normalised energy of the fluctuation structure. In the 
lower panel, the average minor radial (r) locations of low order rational magnetic 
surfaces are shown. 

2.6 An overview of the preprocessed dataset 

As an overview of the preprocessed dataset, figure H] shows the fluctuation 
structures with energy p > 0.2 mapped to the magnetic geometry parameter 
Kh] the radial location of low-order rational magnetic surfaces, t = n/m, are 
also shown. These rational surfaces are important as fluctuations with toroidal 
and poloidal mode numbers n and m respectively can resonate with the twisted 
field lines. The main features of the fluctuation spectra are the resonances 
about Kh = 0.4 and Kh = 0.76, related to the t = 5/4 and t = 4/3 surfaces 
respectively. We expect that any automated process used to locate distinct 
types of fluctuations would identify these features, and hopefully find some 
less obvious features. Indeed in section [Hit can be seen that these two features 
are the first to be distinguished by the following clustering algorithm. 



3 Clustering 

We aim to discover any underlying lower-dimensional model of the dataset; 
that is, groups of fluctuation structures which are similar throughout some 
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range of short time segments. As discussed in section 12.51 we assume that a 
class of fluctuations is localised in the iVc-dimensional Aip-space. For exam- 
ple, it is simple to understand such localisation in terms of a simple cylindri- 
cal geometry with equidistant poloidal measurements, where each mode with 
poloidal mode number m will be located at Atp = 2Trm/Nc in each dimension. 
However, we assume a generalised case in which the fluctuation may have 
arbitrary, including localised, structure. 

Many different types of clustering algorithms exist; here we use the expectation 
maximisation (EM) algorithm which is a method for estimating the most likely 
values of latent variables in a probabilistic model[S]. Here we assume that each 
type of fluctuation can be described by a iVc-dimensional Gaussian distribution 
in Atp space. The latent variables are the mean /i, and standard deviation cTj 
for each cluster i, where i = 1, 2, 3, . . . , Nci and Nci is the number of clusters. 
Given the initial conditions, in the form of random initial /ij and (jj values 
for a prescribed number of clusters, the EM algorithm consists of two steps 
which repeat until a convergence criterion is met. Firstly, the expectation step 
assigns to each datapoint a probability, or expectation value, of belonging 
to each cluster which is calculated with the Gaussian distribution function. 
Secondly, /Xj and cXi are recalculated using the new expectation values as weight 
factors. 

The 10-fold cross-validated log-likelihood ratio is used as a measure of how 
well the cluster assignments fit the data. The cross-validation process involves 
partitioning the dataset into random subsamples and comparing results from 
each subset to avoid oversensitivity to outliers in the data. The likelihood is the 
conditional probability of obtaining the cluster means and standard deviations 
given the observed data. Because the EM algorithm can only guarantee a 
local maximum in likelihood we use a Monte Carlo approach, with multiple 
repetitions with different randomised initial conditions for each Nci- 



3. 1 Visualisation 

The identification of the correct number of clusters, or of those which are 
important, is a task that is by no means trivial to automate. We have found 
inspection of a dendrogram, or cluster tree, mapping to be a practical method 
for identifying the important clusters. The cluster tree displays clusters for 
each Nqi below some maximum value Nci,maxi with all clusters for a given 
Nqi forming a single tree level. Each child cluster is mapped to the cluster on 
the parent level with which it has the largest fraction of common datapoints. 
Cluster branches which do not fork over a significant range of Nci are deemed 
to be well defined, and the point where well defined clusters start to break up 
suggests that Nci is too high. While this procedure is clearly a subjective one. 
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it is effective and does not depend on the type of clustering algorithm used. 

The cluster tree for our example dataset is shown in figures [51 Each cluster 
has been defined only by its phase structure ip and mapped back to Kh and 
frequency / = {2ti)^^uji. The base of the tree {Nci = 1) shows all the data 
within a single cluster (EM: A); as we climb up the tree different classes of 
fluctuation are isolated. For example, the branch starting at cluster EM:B 
contains fluctuations with toroidal mode number n = 5 and poloidal mode 
number m = 4 which occur at configurations near the t = 5/4 resonance 
[nfi ~ 0.4). Similarly, the branch containing cluster EM:C is due to the t = 4/3 
resonance near = 0.75. Other clusters include fluctuations which occur at 
higher order resonances, as well as low frequency n,m = modes (cluster 
EM:0 branch) and weakly defined residual clusters (EM:K and EM:L) which 
would be resolved at a higher level of the tree than is shown here. 

Shown in figure [6] are poloidal phase-angle plots for a single poloidal Mirnov 
array. The centre line corresponds to the cumulative mean phase of the coil 
pairs. The lines above and below are the cumulative cluster standard devia- 
tions of the coil pairs, where (T^ „ = Yl]=i '^jj+i "^i.n = J^jZi Here, 
the magnetic angles have been evaluated for a flux surface at r = 0.1m to 
better represent the broad radial structure expected for these modes. 



4 Discussion 

The data mining algorithm presented here is potentially useful in numerous 
other domains where spatio-temporal data is used. However, there is a limi- 
tation to the nature of the fluctuations amenable to this analysis due to the 
SVD. The SVD is not effective in distinguishing different modes coexisting 
with the same frequency or spatial structure because the modes would share 
a chrono or topo, whereas the SVD requires orthogonal components to dis- 
tinguish modes. The assumption that such coexisting modes are not present 
is also important in assigning a single frequency ui to a fluctuation structure, 
i.e.: two modes with the same spatial structure will also share chronos, but 
only one frequency would be recorded. 

The EM clustering method described here relies on the assumption that clus- 
ters can be described by a Gaussian distribution. To check if the imposed 
Gaussian distributions significantly influence the cluster outcomes, we have 
also used the agglomerative hierarchical (AH) clustering algorithm [0] which 
does not make such an assumption. The initial condition for AH clustering is 
that each fluctuation structure defines a cluster. Using a suitable metric the 
two closest clusters are combined, iterating the process until we have Nci = 1 
gives a (prohibitively large) cluster tree. Compression of the AH cluster tree 
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Fig. 5. Cluster tree of the example dataset. The figure in the bottom left corner is equivalent to figure [H its upper panel shows the 
fluctuation structures mapped to / and k^, the numbers 1 (2000) at the top right are the tree level, Nci , and cluster population 
respectively, EM:A is a cluster label used for reference. For clarity, only a subset of clusters within the tree have their contents displayed 
and EM:G has been displaced to prevent overlap. Vertical parent-child distance is proportional to the distance between cluster means, 
while line thickness is inversely proportional to the Gaussian width of the cluster. Several clusters produced by the agglomerative 
hierarchical (AH: labels) method are also shown, these are essentially equivalent to the Nqi = 10 level EM clusters, see table [1] for 
comparison. 



7r/2 ^ ,'!7r/2 2-k 7r/2 ^ 37r/2 27r 7r/2 ,'!7r/2 27r 

Mean poloidal anqle [r=0.1m] Mean ooloidal anqie [r=0.1m] Mean poloidal anqle [r=0.1m] 

(a) Cluster 47 (b) Cluster 48 (c) cluster 46 

Fig. 6. Variation of phase around one of the poloidal Mirnov arrays, plotted against 
mean magnetic poloidal coordinate. The centre line is the cumulative mean phase of 
the coil pairs, with standard deviation shown above and below. The mode numbers 
shown here are supported by Fourier analysis of the data. 
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A comparison of populations of clusters produced by the EM {Nci = 10) and AH 
algorithms, clusters are shown in figure [5] 



can be achieved by filtering out clusters with small populations, allowing for 
a visualisation similar to the EM cluster tree in figure [5l a set of AH clusters 
which are essentially equivalent to the Nci = 10 level of the EM cluster tree 
are also shown in figure • The clusters resulting from the EM and AH meth- 
ods have been found to be essentially the same apart from the weakly defined 
clusters (EM:K,L) and the 'remainder' (AH:Rem), a quantitative comparison 
between populations of EM and AH clusters in figure [5] is shown in table [TJ 

It is important to consider the scalability and computational requirements 
of the algorithm. Given fixed values of and At, the size of S remains 
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constant and the preprocessing stage has complexity 0{Ns), where Ns is 
the number of timeseries datasets S. The scalabihty of the clustering stage 
depends on the algorithm used, for the EM case we have 0{NciNa), which 
gives 0{Ns) for constant Nci- The AH clustering algorithm is less desirable 
as it has complexity 0{N^) due to distance calculation between each pair of 
fluctuation structures. 

We have implemented the preprocessing and visualisation stages using the 
python language with the Scipy and Matplotlib libraries [7] . The preprocessing 
of our dataset, 4600 S arrays (28 by 1000), takes around 2 hours using a 
1.9 GHz Intel Pentium M processor. The results are stored in MySQL tables; 
a table of fluctuation structure properties excluding Aip-spa.ce mapping is 
around 5 Mb in size, with the 3.6 x 10^ rows of the Aip mapping table taking 
around 30 Mb, using optimal data types. For clustering, we have used the EM 
algorithm from the WEKA suite of data mining tools [8]|9] which runs at about 
0.05x NcixNa CPU seconds using 2.2 GHz AMD Opteron processors. For each 
Nci, 100 randomised initial conditions were used; the results with maximal 
log-likelihood are selected as the best clusters. The WEKA algorithm does not 
operate with toroidal data, so we map the A-^— space from the iVc-dimensional 
torus to a 2A'"c-dimensional cube [—1, 1]2^<= by taking the sm{Ailj) and cos(A^/') 
components. For the present work, no specific efforts were made to optimize 
the clustering process; more efficient clustering routines exist, including, for 
example, genetic algorithms for faster convergence. 

The physical nature of the fluctuations in our example dataset is not yet 
completely understood. The dependence of spectra on plasma density n and 
t suggests a dispersion relation similar to that of the global Alfven eigen- 
mode (GAE) [TUPT] . However, the observed frequencies are smaller than the 
expected GAE frequencies by a factor of around 1/3 [T2|; an experimental 
campaign is presently being undertaken in order to resolve this difference. 



5 Conclusion 

We have presented a highly automated data mining process for the characteri- 
sation of fluctuations in multichannel timeseries data. The manual interaction 
is restricted to two tasks: the selection of a cross-power threshold jmin and 
the choice of appropriate filter parameters. The former requires initial consid- 
eration of threshold effectiveness on a small subset of data, while the latter is 
an operation applied to the dataset whole. 

Given an appropriate choice of clustering algorithm, the data mining process 
scales well, with complexity 0{Ns). We have used the procedure here with 
magnetic fluctuation data from configuration scans in the H-1 heliac, identify- 
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ing different modes in parameter space. The process should be easily adaptable 
to other types of multichannel oscillatory timeseries data. 



Acknowledgements 



The authors would like to thank the H-1 team for continued support of ex- 
perimental operations as well as J. Harris, F. Detering and M. Hegland for 
useful discussions. This work was performed on the H-INF National Plasma 
Fusion Research Facility established by the Australian Government, and oper- 
ated by the Australian National University, with support from the Australian 
Research Council Grant DP0344361 and DP0451960. 



References 



[1] S. M. Hamberger, B. D. Blackwell, L. E. Sharp and D. B. Shenton, H-1 design 
and construction. Fusion Technol. 17 (1990) 123-130 

[2] J. H. Harris et al, Fluctuations and stability of plasmas in the H-lNF heliac. 
Nucl. Fusion 44 (2004) 279-286 

[3] B. D. Blackwell, Results from helical axis stellarators, Phys. Plasmas 8 (2001) 
2238-2244 

[4] T. Dudok de Wit, A.-L. Pecquet, J.-C. Vallet and R. Lima, The biorthogonal 
decomposition as a tool for investigating fluctuations in plasmas. Phys. Plasmas. 
1 (1994) 3288-3300 

[5] A. Dempster, N. Laird and D. Rubin, Maximum likelihood from incomplete 
data via the EM algorithm. J. Royal Stat. Soc. 39 (1977) 1-38 

[6] W. H. E. Day and H. Edelsbrunner, Efficient Algorithms for Agglomerative 
Hierarchical Clustering Methods, J. Classification 1 (1984) 7-24 



[7] |http: / / python.org http://scipy.org, http://matplotlib.sourceforge.net 



[8] I. H. Witten and E. Frank, Data Mining: Practical machine learning tools and 
techniques, 2nd Edition, Morgan Kaufmann (2005) 



[9] http : / / www . cs . waikato. ac . nz/ml/ weka 



[10] K. L. Wong, A review of Alfven eigenmode observations in toroidal plasmas, 
Plasma Phys. Control. Fusion 41 (1999) R1-R56 

[11] D. A. Spong, R. Sanchez and A. Weller, Shear Alfven continua in stellarators, 
Phys. Plasmas 10 (2003) 3217-3224 

[12] D. G. Pretty, PhD Thesis, AustraUan National University (2007) 



14 



